stanEx3
execute mcmc sampling
goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
mcmc=mdl$sample(
data=data,
seed=1,
chains=4,
iter_sampling=smp,
iter_warmup=wrm,
thin=th,
refresh=1000
)
mcmc
}see mcmc result and parameters
seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
print(mcmc)
prs=mcmc$metadata()$model_params[-1] # reject lp__
for(pr in prs){
if(grepl('^y',pr)) next # not show predictive value "y*" information
if(exc!='' && grepl(paste0('^',exc),pr)) next
drw=mcmc$draws(pr)
if(ch){
par(mfrow=c(1,4))
drw[,1,] |> plot(type='l',main='chain1',ylab=pr)
drw[,2,] |> plot(type='l',main='chain2',ylab=pr)
drw[,3,] |> plot(type='l',main='chain3',ylab=pr)
drw[,4,] |> plot(type='l',main='chain4',ylab=pr)
}
par(mfrow=c(1,2))
drw |> hist(main=pr,xlab='')
drw |> density() |> plot(main=pr)
}
}ex8
non linear regression
with one explanatory variable fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)
mcmc=goMCMC(mdl,data,smp,wrm)
mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
}single linear regressionex8-0.stan
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=normal_rng(m1[i],s);
}
}n=20
x=runif(n,0,20)
y=rnorm(n,x*2+5,5)
qplot(x,y)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan')
fn(mdl,data)## Initial log joint probability = -594.909
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 24 -39.8541 0.000404967 0.000456077 1 1 34
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.2 seconds.
## variable estimate
## lp__ -39.85
## b0 4.87
## b1 1.84
## s 4.45
## m0[1] 11.26
## m0[2] 18.67
## m0[3] 9.94
## m0[4] 24.60
## m0[5] 32.34
## m0[6] 24.19
## m0[7] 31.27
## m0[8] 11.26
## m0[9] 9.06
## m0[10] 31.84
## m0[11] 27.11
## m0[12] 30.40
## m0[13] 9.48
## m0[14] 35.05
## m0[15] 36.03
## m0[16] 4.93
## m0[17] 28.48
## m0[18] 33.52
## m0[19] 15.69
## m0[20] 5.45
## y0[1] 15.84
## y0[2] 12.43
## y0[3] 8.42
## y0[4] 22.29
## y0[5] 34.28
## y0[6] 24.69
## y0[7] 31.86
## y0[8] 15.69
## y0[9] 0.22
## y0[10] 30.34
## y0[11] 31.39
## y0[12] 29.04
## y0[13] 2.06
## y0[14] 35.01
## y0[15] 30.44
## y0[16] 6.25
## y0[17] 28.82
## y0[18] 35.81
## y0[19] 19.67
## y0[20] 2.41
## m1[1] 4.93
## m1[2] 8.39
## m1[3] 11.84
## m1[4] 15.30
## m1[5] 18.75
## m1[6] 22.21
## m1[7] 25.66
## m1[8] 29.12
## m1[9] 32.57
## m1[10] 36.03
## y1[1] 3.82
## y1[2] 10.88
## y1[3] 13.01
## y1[4] 5.94
## y1[5] 13.24
## y1[6] 18.10
## y1[7] 27.49
## y1[8] 25.23
## y1[9] 36.65
## y1[10] 31.63
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.4 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -39.95 -39.65 1.23 1.07 -42.36 -38.53 1.00 535 991
## b0 4.96 4.89 2.05 2.15 1.50 8.32 1.00 251 276
## b1 1.83 1.84 0.19 0.19 1.52 2.15 1.00 369 724
## s 5.07 4.93 0.92 0.84 3.78 6.73 1.00 1161 1242
## m0[1] 11.31 11.27 1.55 1.62 8.73 13.90 1.00 265 309
## m0[2] 18.68 18.66 1.19 1.18 16.75 20.67 1.00 429 816
## m0[3] 9.99 9.94 1.64 1.73 7.27 12.71 1.00 259 262
## m0[4] 24.57 24.57 1.21 1.16 22.57 26.52 1.00 1779 1353
## m0[5] 32.27 32.23 1.63 1.56 29.62 35.01 1.00 2559 1522
## m0[6] 24.16 24.16 1.20 1.15 22.21 26.10 1.00 1574 1358
## m0[7] 31.21 31.17 1.56 1.48 28.68 33.81 1.00 2641 1430
## m0[8] 11.31 11.27 1.55 1.62 8.73 13.90 1.00 265 309
## m0[9] 9.12 9.05 1.71 1.80 6.29 11.93 1.00 256 250
## m0[10] 31.77 31.73 1.60 1.53 29.17 34.47 1.00 2597 1498
## m0[11] 27.07 27.05 1.31 1.23 24.95 29.21 1.00 2576 1386
## m0[12] 30.35 30.32 1.50 1.43 27.92 32.83 1.00 2684 1368
## m0[13] 9.54 9.48 1.68 1.77 6.76 12.31 1.00 257 258
## m0[14] 34.96 34.96 1.84 1.75 31.96 38.06 1.00 2204 1471
## m0[15] 35.94 35.94 1.92 1.84 32.76 39.14 1.00 2050 1366
## m0[16] 5.01 4.95 2.04 2.15 1.56 8.36 1.00 251 276
## m0[17] 28.44 28.41 1.38 1.29 26.18 30.71 1.00 2689 1409
## m0[18] 33.44 33.42 1.72 1.65 30.66 36.31 1.00 2442 1451
## m0[19] 15.71 15.68 1.29 1.28 13.61 17.85 1.00 320 588
## m0[20] 5.53 5.47 2.00 2.11 2.16 8.80 1.00 251 267
## y0[1] 11.17 11.03 5.38 4.81 2.19 20.42 1.00 1395 1595
## y0[2] 18.65 18.72 5.28 5.17 10.29 26.93 1.00 1921 1891
## y0[3] 9.96 9.99 5.43 5.33 0.94 18.67 1.00 1322 1839
## y0[4] 24.53 24.53 5.40 5.18 15.55 33.39 1.00 1927 1821
## y0[5] 32.10 32.08 5.42 5.08 23.16 40.95 1.00 1950 1838
## y0[6] 24.16 24.16 5.42 5.14 15.17 32.97 1.00 1936 1787
## y0[7] 31.35 31.42 5.33 5.05 22.65 40.07 1.00 2099 1932
## y0[8] 11.45 11.47 5.37 5.05 2.76 20.16 1.00 1257 1731
## y0[9] 8.97 8.97 5.49 5.18 0.08 18.11 1.00 1448 1573
## y0[10] 31.63 31.69 5.51 5.35 22.71 40.56 1.00 2121 1856
## y0[11] 27.19 27.24 5.42 5.14 18.31 35.90 1.00 1953 1999
## y0[12] 30.31 30.29 5.50 5.18 21.23 39.30 1.00 1955 1913
## y0[13] 9.47 9.41 5.39 5.09 0.68 18.47 1.00 1282 1598
## y0[14] 34.95 34.92 5.29 5.10 26.67 43.65 1.00 2148 1685
## y0[15] 36.05 36.00 5.60 5.21 26.94 45.43 1.00 2045 1855
## y0[16] 5.05 4.86 5.49 5.20 -3.86 14.48 1.00 964 1464
## y0[17] 28.41 28.46 5.31 5.27 19.51 37.15 1.00 2130 1799
## y0[18] 33.22 33.23 5.35 5.13 24.73 41.92 1.00 2271 1889
## y0[19] 15.58 15.60 5.35 5.24 6.78 24.24 1.00 1752 1836
## y0[20] 5.49 5.43 5.62 5.59 -3.84 14.55 1.01 1187 1796
## m1[1] 5.01 4.95 2.04 2.15 1.56 8.36 1.00 251 276
## m1[2] 8.45 8.39 1.76 1.84 5.55 11.37 1.00 254 278
## m1[3] 11.89 11.84 1.51 1.57 9.37 14.43 1.00 268 345
## m1[4] 15.32 15.28 1.31 1.31 13.18 17.49 1.00 312 531
## m1[5] 18.76 18.75 1.19 1.18 16.84 20.74 1.00 434 816
## m1[6] 22.19 22.18 1.16 1.14 20.32 24.11 1.00 835 1059
## m1[7] 25.63 25.62 1.25 1.17 23.59 27.65 1.00 2174 1357
## m1[8] 29.06 29.04 1.42 1.33 26.77 31.39 1.00 2706 1409
## m1[9] 32.50 32.46 1.65 1.58 29.83 35.27 1.00 2539 1523
## m1[10] 35.94 35.94 1.92 1.84 32.76 39.14 1.00 2050 1366
## y1[1] 4.96 5.00 5.52 5.28 -4.41 14.03 1.00 1097 1551
## y1[2] 8.36 8.46 5.32 5.11 -0.51 17.23 1.00 1269 1732
## y1[3] 12.04 12.05 5.11 4.82 4.06 20.46 1.00 1329 1867
## y1[4] 15.27 15.28 5.06 4.84 7.03 23.52 1.00 1971 1894
## y1[5] 18.73 18.83 5.47 5.20 9.61 27.48 1.00 1954 1934
## y1[6] 22.39 22.41 5.26 5.04 13.87 30.75 1.00 1603 1890
## y1[7] 25.59 25.71 5.39 5.29 16.60 34.34 1.00 1950 1891
## y1[8] 29.18 29.33 5.22 4.93 20.54 37.91 1.00 2190 1992
## y1[9] 32.62 32.80 5.31 5.19 23.95 40.86 1.00 2106 1499
## y1[10] 35.98 35.96 5.52 5.26 27.06 45.20 1.00 2011 1941
quadratic regression
y=b0+b2(x-b1)**2 ex8-1.stan
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real b2;
real<lower=0> s;
}
model{
y~normal(b0+b2*(x-b1)^2,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b2*(x[i]-b1)^2;
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b2*(x1[i]-b1)^2;
y1[i]=normal_rng(m1[i],s);
}
}n=20
x=runif(n,0,9)
y=rnorm(n,0.5*(x-4)**2+5,1)
qplot(x,y)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-1.stan')
fn(mdl,data)## Initial log joint probability = -316802
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 99 -35.4836 6.06152 279.657 1 1 199
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 199 -35.4468 0.686777 1461.74 0.2051 0.2051 348
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 299 -35.4432 5.73746 1968.49 1 1 495
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 399 -35.442 0.46783 2597.15 0.3871 0.3871 635
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 499 -35.4412 0.666778 5867.63 0.3691 0.3691 782
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 599 -35.4408 4.5187 4380.86 1 1 934
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 699 -35.4405 0.0200503 262.313 0.153 0.8623 1082
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 799 -35.4403 5.26346 9198.39 0.566 1 1235
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 831 -35.4403 1.52194 625.731 1 1 1276
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -35.44
## b0 166.07
## b1 2585.60
## b2 0.00
## s 3.56
## m0[1] 9.09
## m0[2] 8.67
## m0[3] 8.72
## m0[4] 8.13
## m0[5] 8.80
## m0[6] 8.75
## m0[7] 8.14
## m0[8] 9.06
## m0[9] 8.39
## m0[10] 8.41
## m0[11] 9.07
## m0[12] 9.16
## m0[13] 8.93
## m0[14] 8.22
## m0[15] 8.20
## m0[16] 8.88
## m0[17] 8.82
## m0[18] 8.33
## m0[19] 8.26
## m0[20] 8.47
## y0[1] 10.27
## y0[2] 10.05
## y0[3] 9.86
## y0[4] 0.05
## y0[5] 6.96
## y0[6] 12.64
## y0[7] 10.14
## y0[8] 3.26
## y0[9] 5.32
## y0[10] 6.55
## y0[11] 7.04
## y0[12] 7.29
## y0[13] 7.58
## y0[14] 5.47
## y0[15] 8.17
## y0[16] 9.53
## y0[17] 15.14
## y0[18] 6.79
## y0[19] 11.42
## y0[20] 11.09
## m1[1] 8.13
## m1[2] 8.25
## m1[3] 8.36
## m1[4] 8.47
## m1[5] 8.59
## m1[6] 8.70
## m1[7] 8.82
## m1[8] 8.93
## m1[9] 9.04
## m1[10] 9.16
## y1[1] 15.09
## y1[2] 7.19
## y1[3] 7.32
## y1[4] 12.51
## y1[5] 6.45
## y1[6] 11.89
## y1[7] 9.87
## y1[8] 6.61
## y1[9] 7.75
## y1[10] 8.60
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -12.31 -11.97 1.49 1.23 -15.25 -10.59 1.01 530 461
## b0 4.21 4.20 0.41 0.39 3.55 4.88 1.01 384 498
## b1 3.96 3.96 0.09 0.08 3.82 4.11 1.00 1931 992
## b2 0.56 0.56 0.04 0.04 0.48 0.63 1.01 548 742
## s 1.18 1.15 0.21 0.19 0.88 1.57 1.00 917 866
## m0[1] 13.87 13.90 0.56 0.54 12.92 14.76 1.00 1897 1414
## m0[2] 4.45 4.45 0.41 0.39 3.80 5.12 1.01 391 475
## m0[3] 4.85 4.84 0.40 0.37 4.21 5.50 1.01 407 554
## m0[4] 11.92 11.91 0.59 0.57 10.97 12.86 1.00 1912 1464
## m0[5] 5.99 5.99 0.37 0.33 5.40 6.60 1.01 477 718
## m0[6] 5.26 5.25 0.39 0.35 4.64 5.89 1.01 427 601
## m0[7] 11.70 11.69 0.57 0.56 10.78 12.62 1.00 1952 1513
## m0[8] 12.78 12.80 0.50 0.49 11.94 13.57 1.00 2069 1445
## m0[9] 5.70 5.69 0.35 0.33 5.12 6.28 1.01 514 686
## m0[10] 5.36 5.35 0.36 0.34 4.77 5.95 1.01 466 643
## m0[11] 13.08 13.11 0.52 0.51 12.22 13.91 1.00 2037 1458
## m0[12] 16.51 16.52 0.73 0.70 15.27 17.66 1.00 1276 1250
## m0[13] 8.66 8.67 0.35 0.35 8.09 9.23 1.00 1055 882
## m0[14] 9.28 9.28 0.43 0.41 8.61 9.97 1.00 1949 1403
## m0[15] 9.84 9.83 0.46 0.44 9.11 10.59 1.00 2000 1567
## m0[16] 7.36 7.36 0.35 0.32 6.79 7.92 1.01 656 802
## m0[17] 6.27 6.27 0.36 0.33 5.68 6.86 1.01 503 773
## m0[18] 6.62 6.62 0.35 0.33 6.05 7.18 1.00 728 820
## m0[19] 8.28 8.28 0.38 0.36 7.67 8.91 1.00 1591 1060
## m0[20] 4.68 4.67 0.38 0.36 4.05 5.32 1.01 406 500
## y0[1] 13.87 13.89 1.33 1.26 11.72 16.06 1.00 2109 1654
## y0[2] 4.50 4.48 1.28 1.19 2.38 6.65 1.00 1591 1632
## y0[3] 4.81 4.84 1.28 1.22 2.76 6.91 1.00 1680 1775
## y0[4] 11.89 11.91 1.34 1.30 9.72 14.09 1.00 2011 1711
## y0[5] 5.95 5.94 1.27 1.22 3.85 8.03 1.00 1787 1886
## y0[6] 5.25 5.25 1.24 1.26 3.14 7.29 1.00 1489 1312
## y0[7] 11.75 11.75 1.36 1.30 9.53 13.90 1.00 1546 1593
## y0[8] 12.72 12.71 1.34 1.25 10.48 14.90 1.00 2006 1858
## y0[9] 5.68 5.67 1.29 1.23 3.58 7.78 1.00 1739 1536
## y0[10] 5.37 5.35 1.24 1.22 3.29 7.40 1.00 1874 1819
## y0[11] 13.11 13.10 1.31 1.22 10.97 15.32 1.00 1863 1605
## y0[12] 16.48 16.51 1.36 1.32 14.22 18.61 1.00 1674 1631
## y0[13] 8.66 8.67 1.23 1.20 6.70 10.74 1.00 2028 1878
## y0[14] 9.27 9.25 1.25 1.16 7.21 11.39 1.00 1901 1656
## y0[15] 9.81 9.80 1.27 1.23 7.79 11.89 1.00 2030 1932
## y0[16] 7.36 7.35 1.27 1.23 5.27 9.46 1.00 1861 1860
## y0[17] 6.30 6.30 1.25 1.17 4.26 8.36 1.00 1632 1858
## y0[18] 6.61 6.63 1.24 1.15 4.49 8.61 1.00 1834 1767
## y0[19] 8.31 8.29 1.26 1.22 6.24 10.41 1.00 1815 1839
## y0[20] 4.67 4.63 1.25 1.17 2.68 6.73 1.00 1673 1765
## m1[1] 11.92 11.91 0.59 0.57 10.97 12.86 1.00 1912 1464
## m1[2] 8.53 8.53 0.39 0.38 7.90 9.18 1.00 1710 1179
## m1[3] 6.12 6.11 0.35 0.33 5.55 6.69 1.00 593 759
## m1[4] 4.68 4.67 0.38 0.36 4.05 5.32 1.01 406 500
## m1[5] 4.22 4.21 0.41 0.39 3.56 4.89 1.01 383 501
## m1[6] 4.73 4.72 0.40 0.38 4.08 5.38 1.01 402 528
## m1[7] 6.21 6.21 0.36 0.33 5.62 6.80 1.01 497 768
## m1[8] 8.67 8.67 0.35 0.35 8.09 9.23 1.00 1057 882
## m1[9] 12.10 12.12 0.47 0.46 11.30 12.83 1.00 2115 1425
## m1[10] 16.51 16.52 0.73 0.70 15.27 17.66 1.00 1276 1250
## y1[1] 11.92 11.91 1.36 1.29 9.79 14.18 1.00 2051 1969
## y1[2] 8.55 8.57 1.28 1.21 6.48 10.67 1.00 1875 1854
## y1[3] 6.16 6.18 1.25 1.18 4.06 8.23 1.00 1529 1698
## y1[4] 4.65 4.63 1.21 1.14 2.69 6.63 1.00 1304 1513
## y1[5] 4.15 4.16 1.24 1.23 2.10 6.22 1.00 1417 1736
## y1[6] 4.66 4.68 1.27 1.22 2.56 6.72 1.00 1351 1891
## y1[7] 6.20 6.19 1.28 1.21 4.04 8.27 1.00 1646 1682
## y1[8] 8.65 8.66 1.28 1.29 6.56 10.72 1.00 1724 1728
## y1[9] 12.07 12.10 1.30 1.23 9.94 14.18 1.00 1916 1677
## y1[10] 16.51 16.55 1.41 1.36 14.12 18.77 1.00 1530 1722
both log regression
log y=b0+b1*log x x,y>0
y=exp b0 * x**b1ex8-2.stan
data{
int N;
vector<lower=0>[N] x;
vector<lower=0>[N] y;
int N1;
vector<lower=0>[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~lognormal(b0+b1*log(x),s);
}
generated quantities{
vector[N] m0;
vector<lower=0>[N] y0;
for(i in 1:N){
m0[i]=b0+b1*log(x[i]);
y0[i]=lognormal_rng(m0[i],s);
}
vector[N1] m1;
vector<lower=0>[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*log(x1[i]);
y1[i]=lognormal_rng(m1[i],s);
}
}n=20
x=runif(n,0,10)
y=exp(rnorm(n,log(x)*2+1,1))
grid.arrange(qplot(x,y),
qplot(log(x),log(y)),
ncol=2)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-2.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -67.6061
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 16 -30.3625 0.000408541 0.000581496 1 1 22
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -30.36
## b0 0.78
## b1 2.22
## s 1.10
## m0[1] 5.64
## m0[2] 0.80
## m0[3] 5.07
## m0[4] 4.71
## m0[5] 4.77
## m0[6] 4.63
## m0[7] 5.65
## m0[8] 4.10
## m0[9] 4.29
## m0[10] -0.21
## m0[11] 4.80
## m0[12] 4.77
## m0[13] 5.57
## m0[14] 4.81
## m0[15] 4.52
## m0[16] 3.86
## m0[17] 5.63
## m0[18] 3.20
## m0[19] 2.06
## m0[20] -0.75
## y0[1] 78.58
## y0[2] 0.59
## y0[3] 46.42
## y0[4] 301.26
## y0[5] 122.24
## y0[6] 25.20
## y0[7] 523.62
## y0[8] 114.82
## y0[9] 229.40
## y0[10] 2.21
## y0[11] 61.53
## y0[12] 29.10
## y0[13] 318.82
## y0[14] 161.48
## y0[15] 13.98
## y0[16] 152.85
## y0[17] 164.69
## y0[18] 36.11
## y0[19] 6.96
## y0[20] 0.82
## m1[1] -0.75
## m1[2] 1.59
## m1[3] 2.70
## m1[4] 3.44
## m1[5] 4.00
## m1[6] 4.44
## m1[7] 4.81
## m1[8] 5.13
## m1[9] 5.40
## m1[10] 5.65
## y1[1] 0.56
## y1[2] 40.36
## y1[3] 30.99
## y1[4] 9.95
## y1[5] 80.28
## y1[6] 7.41
## y1[7] 150.57
## y1[8] 48.11
## y1[9] 224.98
## y1[10] 180.53mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.mcmc$metadata()$stan_variables## [1] "lp__" "b0" "b1" "s" "m0" "y0" "m1" "y1"seeMCMC(mcmc,'m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -31.89 -31.54 1.33 1.05 -34.56 -30.47 1.01 772 934
## b0 0.81 0.80 0.56 0.52 -0.08 1.75 1.01 466 545
## b1 2.21 2.21 0.34 0.32 1.63 2.76 1.01 522 729
## s 1.26 1.23 0.23 0.22 0.94 1.68 1.00 1635 1364
## m0[1] 5.64 5.64 0.39 0.37 5.03 6.26 1.00 1684 1289
## m0[2] 0.83 0.82 0.56 0.52 -0.06 1.77 1.01 466 545
## m0[3] 5.07 5.06 0.34 0.32 4.52 5.62 1.00 2109 1384
## m0[4] 4.71 4.71 0.31 0.29 4.21 5.22 1.00 2252 1409
## m0[5] 4.77 4.77 0.31 0.29 4.26 5.29 1.00 2243 1330
## m0[6] 4.63 4.63 0.31 0.28 4.13 5.13 1.00 2253 1430
## m0[7] 5.64 5.64 0.39 0.37 5.03 6.26 1.00 1680 1289
## m0[8] 4.11 4.11 0.29 0.27 3.64 4.57 1.00 1865 1274
## m0[9] 4.29 4.29 0.29 0.27 3.82 4.76 1.00 2110 1381
## m0[10] -0.17 -0.17 0.70 0.66 -1.28 1.03 1.01 458 544
## m0[11] 4.80 4.80 0.32 0.30 4.28 5.33 1.00 2233 1330
## m0[12] 4.77 4.77 0.31 0.29 4.26 5.30 1.00 2242 1330
## m0[13] 5.56 5.56 0.38 0.36 4.96 6.17 1.00 1745 1307
## m0[14] 4.81 4.81 0.32 0.30 4.29 5.34 1.00 2232 1327
## m0[15] 4.52 4.52 0.30 0.28 4.04 5.01 1.00 2236 1409
## m0[16] 3.87 3.87 0.28 0.27 3.40 4.32 1.00 1451 1307
## m0[17] 5.63 5.63 0.39 0.37 5.02 6.25 1.00 1692 1271
## m0[18] 3.21 3.21 0.31 0.29 2.72 3.70 1.01 808 1106
## m0[19] 2.08 2.07 0.40 0.39 1.45 2.77 1.01 520 613
## m0[20] -0.71 -0.72 0.78 0.74 -1.96 0.63 1.01 457 584
## y0[1] 648.83 287.04 1401.59 316.61 28.13 2264.01 1.00 1918 1813
## y0[2] 8.38 2.37 77.23 2.57 0.22 23.92 1.00 1335 1809
## y0[3] 411.92 158.35 1080.41 167.01 18.81 1366.86 1.00 2132 1800
## y0[4] 258.56 111.35 559.60 117.38 12.38 915.14 1.00 2047 1730
## y0[5] 305.95 116.21 1127.08 118.67 14.15 936.63 1.00 1731 1878
## y0[6] 262.61 95.38 776.56 100.43 10.02 817.51 1.00 1909 1806
## y0[7] 743.30 282.54 2518.92 290.81 33.42 2736.99 1.00 1962 1852
## y0[8] 144.97 64.64 280.36 67.75 7.29 530.46 1.00 2006 1837
## y0[9] 198.00 77.98 508.93 83.65 8.62 703.79 1.00 1858 1786
## y0[10] 2.46 0.79 6.74 0.91 0.07 8.91 1.01 1133 1666
## y0[11] 295.55 118.52 833.31 123.40 14.02 965.03 1.00 2074 1816
## y0[12] 344.98 118.19 1333.30 126.62 13.30 1140.70 1.00 1899 1931
## y0[13] 664.06 265.04 1702.03 278.38 28.06 2344.33 1.00 2056 1769
## y0[14] 342.80 127.27 1005.87 134.11 14.15 1068.39 1.00 2010 1655
## y0[15] 251.26 86.80 984.18 90.92 9.76 871.85 1.00 2061 1684
## y0[16] 117.42 48.66 255.17 50.69 5.80 415.75 1.00 2056 2014
## y0[17] 796.30 274.63 3345.23 297.55 32.38 2686.79 1.00 1704 1641
## y0[18] 62.12 24.39 203.79 25.01 2.50 190.41 1.00 2003 1661
## y0[19] 20.85 7.62 48.99 7.96 0.85 77.17 1.00 1850 1931
## y0[20] 1.78 0.48 6.88 0.53 0.04 6.12 1.00 1184 1676
## m1[1] -0.71 -0.72 0.78 0.74 -1.96 0.63 1.01 457 584
## m1[2] 1.61 1.61 0.46 0.43 0.89 2.39 1.01 488 551
## m1[3] 2.72 2.71 0.34 0.33 2.18 3.29 1.01 614 846
## m1[4] 3.45 3.45 0.29 0.28 2.98 3.92 1.00 971 1211
## m1[5] 4.00 4.01 0.29 0.27 3.54 4.46 1.00 1678 1273
## m1[6] 4.44 4.45 0.30 0.27 3.97 4.92 1.00 2203 1388
## m1[7] 4.81 4.81 0.32 0.30 4.29 5.34 1.00 2231 1327
## m1[8] 5.12 5.12 0.34 0.32 4.57 5.68 1.00 2185 1384
## m1[9] 5.40 5.40 0.37 0.34 4.82 5.98 1.00 1903 1287
## m1[10] 5.64 5.64 0.39 0.37 5.03 6.26 1.00 1680 1289
## y1[1] 1.68 0.48 6.20 0.55 0.04 6.36 1.00 1188 1599
## y1[2] 13.45 4.77 69.66 5.28 0.55 41.58 1.00 1570 1724
## y1[3] 40.28 14.95 187.71 15.60 1.72 123.67 1.00 1918 1965
## y1[4] 103.56 32.31 1072.77 34.91 3.54 286.96 1.00 1898 1845
## y1[5] 146.30 54.72 399.81 55.32 7.32 465.89 1.00 2187 1694
## y1[6] 208.31 84.58 486.90 90.75 8.83 737.84 1.00 2036 1947
## y1[7] 313.42 124.30 877.99 132.47 12.99 1136.51 1.00 1899 1855
## y1[8] 421.69 167.36 1030.78 182.31 16.79 1459.20 1.00 2113 1968
## y1[9] 563.86 211.92 1852.93 224.37 24.62 1864.59 1.00 1990 1877
## y1[10] 778.48 279.18 2238.00 299.70 31.75 2664.62 1.00 1892 1766
y0=mcmc$draws('y0')
smy0=summary(y0)
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)
par(mfrow=c(1,2))
hist(log(y)-log(smy0$median),xlab='obs.-prd.',main='residual')
density(log(y)-log(smy0$median)) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=exp(ml5)),xy,col='darkgray')+
geom_line(aes(x=x,y=exp(mu5)),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=exp(m)),xy,col='black')
qplot(log(x),log(y),col=I('red'))+
geom_line(aes(x=log(x),y=ml5),xy,col='darkgray')+
geom_line(aes(x=log(x),y=mu5),xy,col='darkgray')+
geom_line(aes(x=log(x),y=log(yl5)),xy,col='lightgray')+
geom_line(aes(x=log(x),y=log(yu5)),xy,col='lightgray')+
geom_line(aes(x=log(x),y=m),xy,col='black')
ex8-3
exponential increasing/decreasing
y=b0* exp b1* x -> y~N(b0* exp(b1*x),s)
log y=log b0+b1* x -> y~logN(log b0 +b1*x,s)
x,y>0,b0>0
(x=0,y=b0)
b1>0 x->Infinity,y->Infinity
b1<0 x->Infinity,y->+0 n=20
x=runif(n,0,5)
y=rnorm(n,10*exp(-2*x),0.5)
grid.arrange(qplot(x,y),
qplot(x,log(y)),
ncol=2)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)y=b0* exp b1* x -> y~N(b0* exp(b1*x),s) ex8-3-1.stan
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real<lower=0> b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0*exp(b1*x),s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0*exp(b1*x[i]);
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0*exp(b1*x1[i]);
y1[i]=normal_rng(m1[i],s);
}
}mdl=cmdstan_model('./ex8-3-1.stan')
fn(mdl,data)## Initial log joint probability = -94.0922
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 19 9.308 7.50058e-05 0.000441256 1 1 40
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ 9.31
## b0 8.95
## b1 -2.06
## s 0.38
## m0[1] 0.01
## m0[2] 0.01
## m0[3] 0.01
## m0[4] 0.01
## m0[5] 0.00
## m0[6] 3.26
## m0[7] 1.63
## m0[8] 0.04
## m0[9] 0.00
## m0[10] 0.00
## m0[11] 1.20
## m0[12] 3.93
## m0[13] 0.00
## m0[14] 0.01
## m0[15] 2.13
## m0[16] 1.68
## m0[17] 0.00
## m0[18] 1.43
## m0[19] 0.39
## m0[20] 0.00
## y0[1] 0.04
## y0[2] -0.83
## y0[3] 0.01
## y0[4] 0.63
## y0[5] 0.53
## y0[6] 3.40
## y0[7] 1.38
## y0[8] -0.68
## y0[9] 0.29
## y0[10] -0.45
## y0[11] 1.52
## y0[12] 4.11
## y0[13] 0.03
## y0[14] -0.69
## y0[15] 1.56
## y0[16] 1.63
## y0[17] -0.35
## y0[18] 1.77
## y0[19] 0.34
## y0[20] 0.15
## m1[1] 3.93
## m1[2] 1.46
## m1[3] 0.54
## m1[4] 0.20
## m1[5] 0.07
## m1[6] 0.03
## m1[7] 0.01
## m1[8] 0.00
## m1[9] 0.00
## m1[10] 0.00
## y1[1] 4.18
## y1[2] 1.58
## y1[3] 0.19
## y1[4] 0.05
## y1[5] 0.79
## y1[6] -0.32
## y1[7] 0.29
## y1[8] -0.39
## y1[9] 0.03
## y1[10] -0.40
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ 10.51 10.83 1.31 1.14 7.94 11.95 1.01 634 781
## b0 9.66 9.41 2.14 1.65 6.93 13.44 1.00 704 716
## b1 -2.18 -2.14 0.35 0.31 -2.79 -1.66 1.00 751 874
## s 0.43 0.42 0.08 0.08 0.32 0.58 1.01 801 867
## m0[1] 0.01 0.01 0.01 0.01 0.00 0.03 1.00 780 870
## m0[2] 0.01 0.01 0.01 0.00 0.00 0.02 1.00 778 885
## m0[3] 0.02 0.01 0.01 0.01 0.00 0.04 1.00 781 870
## m0[4] 0.01 0.00 0.01 0.00 0.00 0.02 1.00 778 885
## m0[5] 0.00 0.00 0.00 0.00 0.00 0.01 1.00 772 844
## m0[6] 3.27 3.26 0.24 0.23 2.87 3.67 1.00 1363 1181
## m0[7] 1.58 1.58 0.19 0.18 1.27 1.88 1.00 1215 1123
## m0[8] 0.04 0.03 0.03 0.02 0.01 0.09 1.00 790 891
## m0[9] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 770 843
## m0[10] 0.00 0.00 0.00 0.00 0.00 0.01 1.00 774 885
## m0[11] 1.15 1.14 0.18 0.17 0.85 1.44 1.00 1038 998
## m0[12] 3.98 3.96 0.35 0.32 3.44 4.58 1.00 986 1137
## m0[13] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 771 843
## m0[14] 0.01 0.01 0.01 0.01 0.00 0.04 1.00 781 870
## m0[15] 2.09 2.08 0.19 0.17 1.78 2.39 1.00 1546 1112
## m0[16] 1.63 1.63 0.19 0.18 1.32 1.93 1.00 1243 1122
## m0[17] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 771 843
## m0[18] 1.38 1.38 0.19 0.18 1.07 1.68 1.00 1122 1058
## m0[19] 0.37 0.36 0.12 0.11 0.19 0.57 1.00 865 815
## m0[20] 0.00 0.00 0.00 0.00 0.00 0.01 1.00 774 885
## y0[1] 0.02 0.03 0.45 0.44 -0.71 0.75 1.00 1832 1777
## y0[2] 0.01 0.01 0.44 0.42 -0.71 0.73 1.00 1928 1799
## y0[3] 0.03 0.04 0.44 0.43 -0.69 0.74 1.00 1946 1675
## y0[4] 0.01 0.01 0.44 0.43 -0.71 0.76 1.00 1685 1831
## y0[5] -0.02 -0.02 0.44 0.41 -0.75 0.70 1.00 1960 1776
## y0[6] 3.28 3.28 0.51 0.49 2.47 4.10 1.00 1646 1506
## y0[7] 1.59 1.60 0.48 0.46 0.79 2.36 1.00 1935 1847
## y0[8] 0.05 0.06 0.44 0.44 -0.67 0.74 1.00 1914 1673
## y0[9] 0.01 0.01 0.45 0.42 -0.70 0.77 1.00 2005 1872
## y0[10] 0.00 0.01 0.44 0.42 -0.72 0.71 1.00 1866 1703
## y0[11] 1.15 1.15 0.49 0.46 0.35 1.94 1.00 1815 1771
## y0[12] 3.96 3.96 0.56 0.55 3.10 4.89 1.00 1399 1178
## y0[13] -0.01 -0.01 0.43 0.40 -0.71 0.68 1.00 1791 1870
## y0[14] 0.02 0.03 0.43 0.40 -0.68 0.73 1.00 1920 1768
## y0[15] 2.07 2.09 0.47 0.43 1.26 2.84 1.00 1927 1730
## y0[16] 1.61 1.61 0.48 0.44 0.83 2.40 1.00 1800 1846
## y0[17] 0.00 0.00 0.43 0.42 -0.73 0.69 1.00 2048 1894
## y0[18] 1.38 1.37 0.47 0.43 0.62 2.18 1.00 1924 1723
## y0[19] 0.36 0.36 0.46 0.43 -0.38 1.11 1.00 1835 1684
## y0[20] -0.01 -0.01 0.44 0.42 -0.74 0.69 1.00 2024 1814
## m1[1] 3.98 3.96 0.35 0.32 3.44 4.58 1.00 986 1137
## m1[2] 1.40 1.40 0.19 0.18 1.10 1.70 1.00 1133 1058
## m1[3] 0.51 0.50 0.14 0.13 0.29 0.75 1.00 889 882
## m1[4] 0.19 0.18 0.08 0.07 0.08 0.33 1.00 824 811
## m1[5] 0.07 0.06 0.04 0.04 0.02 0.15 1.00 798 891
## m1[6] 0.03 0.02 0.02 0.02 0.01 0.07 1.00 787 902
## m1[7] 0.01 0.01 0.01 0.01 0.00 0.03 1.00 780 870
## m1[8] 0.00 0.00 0.01 0.00 0.00 0.01 1.00 775 885
## m1[9] 0.00 0.00 0.00 0.00 0.00 0.01 1.00 772 885
## m1[10] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 771 843
## y1[1] 3.98 3.96 0.57 0.54 3.07 4.93 1.00 1374 1467
## y1[2] 1.41 1.40 0.48 0.46 0.61 2.22 1.00 1661 1645
## y1[3] 0.50 0.50 0.46 0.44 -0.26 1.25 1.00 1784 1861
## y1[4] 0.19 0.18 0.45 0.42 -0.53 0.92 1.00 1954 1755
## y1[5] 0.08 0.08 0.44 0.43 -0.64 0.83 1.00 2001 1771
## y1[6] 0.03 0.03 0.44 0.41 -0.68 0.76 1.00 2223 1849
## y1[7] 0.00 -0.01 0.43 0.41 -0.69 0.71 1.00 1926 1591
## y1[8] 0.00 0.01 0.42 0.40 -0.71 0.68 1.00 2282 1942
## y1[9] -0.02 -0.02 0.43 0.41 -0.75 0.70 1.00 2001 1638
## y1[10] 0.00 0.01 0.44 0.43 -0.75 0.74 1.00 1928 1838
log y=log b0+b1* x -> y~logN(log b0 +b1*x,s)ex8-3-2.stan
data{
int N;
vector<lower=0>[N] x;
vector<lower=0>[N] y;
int N1;
vector[N1] x1;
}
parameters{
real<lower=0> b0;
real<lower=-10,upper=10> b1;
real<lower=0> s;
}
model{
y~lognormal(log(b0)+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=log(b0)+b1*x[i];
y0[i]=lognormal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=log(b0)+b1*x1[i];
y1[i]=lognormal_rng(m1[i],s);
}
}n=20
x=runif(n,0,5)
y=rlnorm(n,log(10)-2*x,0.5)
grid.arrange(qplot(x,y),
qplot(x,log(y)),
ncol=2)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-3-2.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -17976.8
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 48 -9.81683 1.26169e-05 0.000321397 0.2257 0.8647 135
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -9.82
## b0 9.53
## b1 -1.97
## s 0.40
## m0[1] 1.72
## m0[2] -1.69
## m0[3] -4.29
## m0[4] 2.02
## m0[5] 1.89
## m0[6] -2.69
## m0[7] -5.32
## m0[8] -1.51
## m0[9] -4.96
## m0[10] -0.48
## m0[11] -5.82
## m0[12] -0.52
## m0[13] -5.40
## m0[14] 1.70
## m0[15] 1.38
## m0[16] 0.81
## m0[17] -1.98
## m0[18] -6.40
## m0[19] -4.15
## m0[20] 1.49
## y0[1] 7.86
## y0[2] 0.09
## y0[3] 0.01
## y0[4] 4.56
## y0[5] 9.98
## y0[6] 0.03
## y0[7] 0.01
## y0[8] 0.17
## y0[9] 0.01
## y0[10] 1.05
## y0[11] 0.00
## y0[12] 1.02
## y0[13] 0.00
## y0[14] 6.44
## y0[15] 6.24
## y0[16] 3.84
## y0[17] 0.21
## y0[18] 0.00
## y0[19] 0.01
## y0[20] 4.81
## m1[1] 2.02
## m1[2] 1.08
## m1[3] 0.15
## m1[4] -0.79
## m1[5] -1.72
## m1[6] -2.66
## m1[7] -3.59
## m1[8] -4.53
## m1[9] -5.46
## m1[10] -6.40
## y1[1] 7.11
## y1[2] 3.17
## y1[3] 0.59
## y1[4] 1.15
## y1[5] 0.15
## y1[6] 0.08
## y1[7] 0.02
## y1[8] 0.01
## y1[9] 0.00
## y1[10] 0.00mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.mcmc$metadata()$stan_variables## [1] "lp__" "b0" "b1" "s" "m0" "y0" "m1" "y1"seeMCMC(mcmc,'m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -8.53 -8.16 1.38 1.05 -10.95 -7.09 1.01 433 545
## b0 9.96 9.75 1.86 1.75 7.22 13.27 1.01 597 446
## b1 -1.98 -1.98 0.07 0.07 -2.10 -1.87 1.01 748 720
## s 0.45 0.44 0.08 0.07 0.34 0.60 1.01 698 609
## m0[1] 1.75 1.74 0.17 0.16 1.47 2.03 1.01 596 425
## m0[2] -1.68 -1.68 0.10 0.10 -1.85 -1.51 1.00 1078 1363
## m0[3] -4.29 -4.29 0.13 0.13 -4.50 -4.07 1.01 2155 1193
## m0[4] 2.04 2.04 0.18 0.17 1.75 2.34 1.01 597 437
## m0[5] 1.91 1.91 0.18 0.17 1.63 2.20 1.01 597 431
## m0[6] -2.68 -2.68 0.11 0.10 -2.86 -2.51 1.00 1777 1166
## m0[7] -5.33 -5.33 0.16 0.15 -5.59 -5.07 1.00 1917 1213
## m0[8] -1.50 -1.50 0.11 0.10 -1.67 -1.33 1.00 990 1303
## m0[9] -4.96 -4.97 0.15 0.14 -5.21 -4.72 1.01 2010 1159
## m0[10] -0.47 -0.46 0.12 0.11 -0.66 -0.28 1.01 709 747
## m0[11] -5.83 -5.84 0.17 0.17 -6.12 -5.55 1.00 1764 1238
## m0[12] -0.50 -0.50 0.12 0.11 -0.70 -0.32 1.01 714 747
## m0[13] -5.41 -5.41 0.16 0.16 -5.67 -5.15 1.00 1897 1197
## m0[14] 1.72 1.72 0.17 0.16 1.45 2.01 1.01 597 425
## m0[15] 1.40 1.40 0.16 0.16 1.14 1.68 1.01 601 460
## m0[16] 0.83 0.83 0.15 0.14 0.59 1.07 1.01 613 513
## m0[17] -1.97 -1.97 0.10 0.10 -2.15 -1.81 1.00 1259 1347
## m0[18] -6.41 -6.42 0.19 0.19 -6.72 -6.10 1.00 1593 1167
## m0[19] -4.16 -4.16 0.13 0.13 -4.36 -3.94 1.01 2172 1109
## m0[20] 1.52 1.51 0.17 0.16 1.25 1.79 1.01 599 459
## y0[1] 6.54 5.90 3.51 2.67 2.58 12.76 1.00 1532 1513
## y0[2] 0.21 0.19 0.10 0.08 0.09 0.40 1.00 1993 1645
## y0[3] 0.02 0.01 0.01 0.01 0.01 0.03 1.00 2013 2023
## y0[4] 8.72 7.70 4.78 3.58 3.43 16.92 1.00 1931 1727
## y0[5] 7.50 6.68 3.88 3.15 3.01 14.60 1.00 1528 1688
## y0[6] 0.08 0.07 0.04 0.03 0.03 0.15 1.00 2129 1840
## y0[7] 0.01 0.00 0.00 0.00 0.00 0.01 1.00 1877 1845
## y0[8] 0.25 0.22 0.12 0.10 0.10 0.46 1.00 1843 1739
## y0[9] 0.01 0.01 0.00 0.00 0.00 0.01 1.00 1684 1741
## y0[10] 0.71 0.63 0.38 0.28 0.29 1.39 1.00 2158 1994
## y0[11] 0.00 0.00 0.00 0.00 0.00 0.01 1.00 1939 1970
## y0[12] 0.66 0.59 0.33 0.26 0.28 1.27 1.00 1824 2015
## y0[13] 0.01 0.00 0.00 0.00 0.00 0.01 1.00 2110 2042
## y0[14] 6.35 5.64 3.30 2.57 2.63 12.67 1.00 1662 1690
## y0[15] 4.67 4.07 2.48 1.85 1.90 9.33 1.00 1438 1332
## y0[16] 2.52 2.25 1.28 1.01 1.01 4.99 1.00 1642 1443
## y0[17] 0.16 0.14 0.08 0.06 0.06 0.31 1.00 1983 1940
## y0[18] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 1939 1858
## y0[19] 0.02 0.02 0.01 0.01 0.01 0.04 1.00 2195 1922
## y0[20] 5.06 4.42 2.63 2.01 2.04 9.72 1.00 1545 1584
## m1[1] 2.04 2.04 0.18 0.17 1.75 2.34 1.01 597 437
## m1[2] 1.10 1.10 0.15 0.15 0.86 1.36 1.01 606 498
## m1[3] 0.16 0.16 0.13 0.13 -0.05 0.37 1.01 647 584
## m1[4] -0.78 -0.77 0.11 0.11 -0.96 -0.60 1.00 762 870
## m1[5] -1.72 -1.72 0.10 0.10 -1.89 -1.55 1.00 1097 1379
## m1[6] -2.66 -2.65 0.11 0.10 -2.83 -2.48 1.00 1761 1166
## m1[7] -3.60 -3.59 0.12 0.12 -3.79 -3.41 1.01 2158 1211
## m1[8] -4.54 -4.54 0.14 0.14 -4.76 -4.31 1.01 2113 1193
## m1[9] -5.47 -5.48 0.16 0.16 -5.74 -5.21 1.00 1877 1197
## m1[10] -6.41 -6.42 0.19 0.19 -6.72 -6.10 1.00 1593 1167
## y1[1] 8.74 7.89 4.28 3.53 3.54 16.98 1.00 1838 1896
## y1[2] 3.39 3.06 1.74 1.31 1.38 6.38 1.00 1836 1512
## y1[3] 1.34 1.21 0.69 0.54 0.54 2.54 1.00 1546 1816
## y1[4] 0.50 0.45 0.26 0.19 0.22 0.95 1.00 1972 1884
## y1[5] 0.21 0.18 0.11 0.08 0.09 0.41 1.00 1827 1920
## y1[6] 0.08 0.07 0.04 0.03 0.03 0.15 1.00 1931 1804
## y1[7] 0.03 0.03 0.02 0.01 0.01 0.06 1.00 1983 2014
## y1[8] 0.01 0.01 0.01 0.00 0.01 0.02 1.00 2046 1753
## y1[9] 0.00 0.00 0.00 0.00 0.00 0.01 1.00 1707 1709
## y1[10] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 2272 1860
y0=mcmc$draws('y0')
smy0=summary(y0)
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,log(y),col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=log(yl5)),xy,col='lightgray')+
geom_line(aes(x=x,y=log(yu5)),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
ex8-4
growth curve
y=b0* (1-exp(-b1* x)) -> y~N(1-exp(-b1*x),s)
x,y>0, b0,b1>0
(x=0,y=0), (x->Infinity,y->b0)data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real<lower=0,upper=100> b0;
real<lower=0,upper=10> b1;
real<lower=0> s;
}
model{
y~normal(b0*(1-exp(-b1*x)),s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0*(1-exp(-b1*x[i]));
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0*(1-exp(-b1*x1[i]));
y1[i]=normal_rng(m1[i],s);
}
}n=20
x=runif(n,0,9)
y=rnorm(n,10*(1-exp(-0.5*x)),1)
qplot(x,y)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-4.stan')
fn(mdl,data)## Initial log joint probability = -1866.79
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 26 -9.36325 0.000146461 0.00351921 0.9112 0.9112 60
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -9.36
## b0 11.49
## b1 0.36
## s 0.97
## m0[1] 8.15
## m0[2] 3.43
## m0[3] 9.88
## m0[4] 10.42
## m0[5] 1.17
## m0[6] 10.72
## m0[7] 9.47
## m0[8] 2.22
## m0[9] 5.98
## m0[10] 11.02
## m0[11] 8.43
## m0[12] 10.27
## m0[13] 10.03
## m0[14] 8.28
## m0[15] 2.28
## m0[16] 10.11
## m0[17] 10.98
## m0[18] 1.17
## m0[19] 9.36
## m0[20] 11.00
## y0[1] 9.89
## y0[2] 1.47
## y0[3] 9.72
## y0[4] 10.42
## y0[5] 0.81
## y0[6] 11.13
## y0[7] 9.08
## y0[8] 1.46
## y0[9] 6.07
## y0[10] 11.64
## y0[11] 8.58
## y0[12] 10.48
## y0[13] 10.74
## y0[14] 7.42
## y0[15] 2.36
## y0[16] 7.12
## y0[17] 10.54
## y0[18] 0.84
## y0[19] 10.56
## y0[20] 10.19
## m1[1] 1.17
## m1[2] 4.17
## m1[3] 6.30
## m1[4] 7.81
## m1[5] 8.88
## m1[6] 9.64
## m1[7] 10.18
## m1[8] 10.56
## m1[9] 10.83
## m1[10] 11.02
## y1[1] 2.04
## y1[2] 4.21
## y1[3] 6.93
## y1[4] 6.87
## y1[5] 6.87
## y1[6] 9.52
## y1[7] 11.19
## y1[8] 9.10
## y1[9] 11.35
## y1[10] 10.78
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -9.76 -9.40 1.36 1.11 -12.55 -8.29 1.00 567 552
## b0 11.51 11.43 0.81 0.72 10.36 12.95 1.00 678 677
## b1 0.38 0.37 0.08 0.07 0.27 0.53 1.00 598 617
## s 1.10 1.07 0.20 0.19 0.82 1.47 1.00 871 747
## m0[1] 8.18 8.19 0.42 0.41 7.49 8.84 1.00 890 1162
## m0[2] 3.50 3.47 0.41 0.37 2.89 4.25 1.00 655 537
## m0[3] 9.86 9.87 0.30 0.28 9.38 10.32 1.00 2104 1269
## m0[4] 10.38 10.39 0.33 0.32 9.84 10.93 1.00 1382 1083
## m0[5] 1.21 1.19 0.18 0.16 0.95 1.53 1.00 639 494
## m0[6] 10.68 10.69 0.39 0.37 10.03 11.32 1.00 992 889
## m0[7] 9.47 9.47 0.31 0.29 8.93 9.95 1.00 1822 1286
## m0[8] 2.28 2.25 0.30 0.27 1.84 2.84 1.00 646 494
## m0[9] 6.05 6.03 0.51 0.47 5.25 6.93 1.00 704 660
## m0[10] 10.98 10.99 0.49 0.46 10.19 11.79 1.00 810 764
## m0[11] 8.45 8.46 0.40 0.38 7.79 9.07 1.00 962 1151
## m0[12] 10.24 10.24 0.31 0.30 9.73 10.75 1.00 1671 1028
## m0[13] 10.00 10.01 0.30 0.28 9.52 10.48 1.00 2058 1102
## m0[14] 8.30 8.31 0.41 0.39 7.63 8.95 1.00 920 1191
## m0[15] 2.34 2.31 0.31 0.28 1.89 2.91 1.00 647 494
## m0[16] 10.09 10.09 0.30 0.28 9.59 10.57 1.00 2000 1042
## m0[17] 10.94 10.95 0.47 0.45 10.17 11.73 1.00 828 847
## m0[18] 1.21 1.19 0.17 0.16 0.95 1.53 1.00 639 494
## m0[19] 9.36 9.36 0.32 0.30 8.80 9.84 1.00 1650 1323
## m0[20] 10.96 10.97 0.48 0.46 10.18 11.76 1.00 818 764
## y0[1] 8.20 8.19 1.21 1.18 6.24 10.17 1.00 1935 1968
## y0[2] 3.55 3.53 1.23 1.18 1.67 5.69 1.00 1642 1913
## y0[3] 9.86 9.86 1.17 1.10 7.99 11.78 1.00 1910 1833
## y0[4] 10.36 10.35 1.14 1.08 8.51 12.25 1.00 1937 1881
## y0[5] 1.20 1.21 1.12 1.06 -0.67 3.01 1.00 1930 1915
## y0[6] 10.68 10.71 1.22 1.12 8.67 12.58 1.00 1899 1786
## y0[7] 9.47 9.48 1.14 1.06 7.63 11.30 1.00 1617 1805
## y0[8] 2.29 2.30 1.18 1.11 0.38 4.23 1.00 1799 1394
## y0[9] 6.07 6.05 1.24 1.19 4.08 8.17 1.00 1583 1750
## y0[10] 10.98 11.01 1.25 1.17 8.94 13.05 1.00 1504 1824
## y0[11] 8.43 8.45 1.17 1.17 6.55 10.32 1.00 1727 1958
## y0[12] 10.22 10.22 1.16 1.10 8.27 12.09 1.00 1785 1641
## y0[13] 10.01 10.00 1.16 1.12 8.12 11.99 1.00 2165 1954
## y0[14] 8.29 8.28 1.16 1.09 6.37 10.17 1.00 1838 1856
## y0[15] 2.33 2.28 1.15 1.14 0.52 4.33 1.00 1686 1891
## y0[16] 10.07 10.07 1.16 1.13 8.11 11.94 1.00 1910 1956
## y0[17] 10.95 11.00 1.23 1.21 8.90 12.92 1.00 1938 1894
## y0[18] 1.21 1.20 1.16 1.08 -0.71 3.07 1.00 1886 1820
## y0[19] 9.32 9.32 1.17 1.10 7.38 11.30 1.00 1980 1933
## y0[20] 10.94 10.93 1.23 1.21 8.96 12.94 1.00 1611 1514
## m1[1] 1.21 1.19 0.17 0.16 0.95 1.53 1.00 639 494
## m1[2] 4.25 4.22 0.46 0.42 3.56 5.07 1.00 661 527
## m1[3] 6.37 6.35 0.50 0.47 5.57 7.23 1.00 714 725
## m1[4] 7.85 7.84 0.44 0.43 7.12 8.56 1.00 826 951
## m1[5] 8.89 8.90 0.36 0.33 8.27 9.45 1.00 1162 1341
## m1[6] 9.62 9.63 0.30 0.28 9.10 10.10 1.00 1993 1323
## m1[7] 10.15 10.15 0.31 0.29 9.65 10.64 1.00 1865 1164
## m1[8] 10.52 10.52 0.36 0.34 9.94 11.09 1.00 1171 928
## m1[9] 10.79 10.80 0.42 0.40 10.11 11.48 1.00 908 859
## m1[10] 10.98 10.99 0.49 0.46 10.19 11.79 1.00 810 764
## y1[1] 1.23 1.22 1.14 1.12 -0.62 3.13 1.00 1877 1842
## y1[2] 4.25 4.27 1.19 1.11 2.34 6.21 1.00 1748 1522
## y1[3] 6.36 6.35 1.26 1.19 4.34 8.41 1.00 1601 1772
## y1[4] 7.86 7.85 1.22 1.17 5.87 9.82 1.00 1853 1575
## y1[5] 8.89 8.90 1.14 1.06 7.00 10.68 1.00 2045 1900
## y1[6] 9.63 9.66 1.17 1.16 7.77 11.49 1.00 1930 1772
## y1[7] 10.18 10.18 1.16 1.09 8.33 12.02 1.00 2176 1897
## y1[8] 10.51 10.54 1.19 1.09 8.55 12.42 1.00 1801 1856
## y1[9] 10.80 10.79 1.16 1.12 8.93 12.71 1.00 1849 991
## y1[10] 11.01 11.00 1.19 1.15 9.08 12.99 1.00 1692 2012
ex8-5
sigmoid curve
y=Ym/ 1+exp(-b1* (x-b0)) -> y~B(Ym, 1+exp(-b1*(x-b0)))
b0,b1>0
x[0,Infinity), y[0,Ym]
(x=b0, y=Ym/2)data{
int N;
vector[N] x;
int Ym;
array[N] int y;
int N1;
vector[N1] x1;
}
parameters{
real<lower=0,upper=100> b0;
real<lower=0,upper=100> b1;
}
model{
y~binomial_logit(Ym,b1*(x-b0));
}
generated quantities{
array[N] int y0;
for(i in 1:N){
y0[i]=binomial_rng(Ym,inv_logit(b1*x[i]-b0*b1));
}
array[N1] int y1;
for(i in 1:N1){
y1[i]=binomial_rng(Ym,inv_logit(b1*x1[i]-b0*b1));
}
}n=20
x=runif(n,0,9)
ym=10
y=rbinom(n,ym,1/(1+exp(-2*(x-4))))
qplot(x,y)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,Ym=ym,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-5.stan')
mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.mcmc$metadata()$stan_variables## [1] "lp__" "b0" "b1" "y0" "y1"seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -16.61 -16.31 1.03 0.72 -18.66 -15.65 1.01 626 584
## b0 4.00 3.98 0.26 0.26 3.59 4.44 1.00 1191 1273
## b1 2.08 2.01 0.47 0.43 1.45 2.92 1.01 378 257
## y0[1] 0.09 0.00 0.31 0.00 0.00 1.00 1.00 2023 2027
## y0[2] 2.68 3.00 1.61 1.48 0.00 6.00 1.00 1884 2041
## y0[3] 10.00 10.00 0.06 0.00 10.00 10.00 1.00 2030 NA
## y0[4] 10.00 10.00 0.07 0.00 10.00 10.00 1.00 1700 NA
## y0[5] 9.97 10.00 0.16 0.00 10.00 10.00 1.00 1681 NA
## y0[6] 0.31 0.00 0.61 0.00 0.00 2.00 1.00 1500 1462
## y0[7] 0.25 0.00 0.53 0.00 0.00 1.00 1.00 1797 1824
## y0[8] 9.93 10.00 0.28 0.00 9.00 10.00 1.00 1589 NA
## y0[9] 0.24 0.00 0.50 0.00 0.00 1.00 1.00 1707 1754
## y0[10] 9.83 10.00 0.43 0.00 9.00 10.00 1.00 1491 NA
## y0[11] 10.00 10.00 0.06 0.00 10.00 10.00 1.00 2026 NA
## y0[12] 9.96 10.00 0.20 0.00 10.00 10.00 1.00 1607 NA
## y0[13] 9.96 10.00 0.20 0.00 10.00 10.00 1.00 1922 NA
## y0[14] 9.84 10.00 0.44 0.00 9.00 10.00 1.00 955 NA
## y0[15] 1.87 2.00 1.41 1.48 0.00 4.00 1.00 1913 1915
## y0[16] 10.00 10.00 0.04 0.00 10.00 10.00 1.00 2023 NA
## y0[17] 0.07 0.00 0.27 0.00 0.00 1.00 1.00 1576 1599
## y0[18] 9.90 10.00 0.33 0.00 9.00 10.00 1.00 1440 NA
## y0[19] 10.00 10.00 0.06 0.00 10.00 10.00 1.00 1655 NA
## y0[20] 0.01 0.00 0.11 0.00 0.00 0.00 1.00 1744 1744
## y1[1] 0.01 0.00 0.08 0.00 0.00 0.00 1.00 2043 2043
## y1[2] 0.05 0.00 0.22 0.00 0.00 0.00 1.00 1747 1777
## y1[3] 0.26 0.00 0.53 0.00 0.00 1.00 1.00 1739 1711
## y1[4] 1.43 1.00 1.21 1.48 0.00 4.00 1.00 1976 1785
## y1[5] 5.18 5.00 1.94 1.48 2.00 8.00 1.00 1548 1548
## y1[6] 8.59 9.00 1.38 1.48 6.00 10.00 1.00 1011 NA
## y1[7] 9.69 10.00 0.63 0.00 8.00 10.00 1.00 1257 NA
## y1[8] 9.93 10.00 0.27 0.00 9.00 10.00 1.00 1600 NA
## y1[9] 9.98 10.00 0.13 0.00 10.00 10.00 1.00 1861 NA
## y1[10] 10.00 10.00 0.06 0.00 10.00 10.00 1.00 2029 NA
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,ymed=smy1$median,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=ymed),xy,col='black')
ex8-6
up and down
y=b0(exp(-b1* x)-exp(-b2* x)) -> y~N(b0(exp(-b1* x)-exp(-b2* x)),s)
b0,b1,b2>0, b1<b2
x[0,Infinity), 0<y<b0
(x=log b1-log b2 / b1-b2, y=max(y))data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real<lower=0,upper=200> b0;
real<lower=0,upper=1> b1;
real<lower=0,upper=1> b2;
real<lower=0,upper=10> s;
}
model{
y~normal(b0*(exp(-b1*x)-exp(-b2*x)),s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0*(exp(-b1*x[i])-exp(-b2*x[i]));
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0*(exp(-b1*x1[i])-exp(-b2*x1[i]));
y1[i]=normal_rng(m1[i],s);
}
}n=20
x=runif(n,0,50)
y=rnorm(n,100*(exp(-0.03*x)-exp(-0.2*x)),1)
qplot(x,y)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-6.stan')
fn(mdl,data)## Initial log joint probability = -1235.08
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 69 3.80917 1.97211e-05 0.0160704 1 1 117
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ 3.81
## b0 104.77
## b1 0.03
## b2 0.19
## s 0.50
## m0[1] 57.76
## m0[2] 60.85
## m0[3] 38.52
## m0[4] 60.80
## m0[5] 32.73
## m0[6] 46.81
## m0[7] 61.40
## m0[8] 61.41
## m0[9] 22.94
## m0[10] 28.87
## m0[11] 37.08
## m0[12] 39.17
## m0[13] 24.72
## m0[14] 49.48
## m0[15] 23.99
## m0[16] 57.04
## m0[17] 54.01
## m0[18] 21.98
## m0[19] 44.47
## m0[20] 60.73
## y0[1] 57.69
## y0[2] 60.67
## y0[3] 38.61
## y0[4] 61.40
## y0[5] 32.52
## y0[6] 46.53
## y0[7] 61.76
## y0[8] 61.27
## y0[9] 22.87
## y0[10] 28.73
## y0[11] 37.31
## y0[12] 38.96
## y0[13] 23.74
## y0[14] 49.43
## y0[15] 24.21
## y0[16] 56.89
## y0[17] 54.50
## y0[18] 21.89
## y0[19] 44.42
## y0[20] 60.66
## m1[1] 28.87
## m1[2] 57.76
## m1[3] 61.06
## m1[4] 56.14
## m1[5] 49.17
## m1[6] 42.24
## m1[7] 36.00
## m1[8] 30.57
## m1[9] 25.93
## m1[10] 21.98
## y1[1] 29.51
## y1[2] 57.63
## y1[3] 61.11
## y1[4] 56.19
## y1[5] 48.70
## y1[6] 42.11
## y1[7] 36.03
## y1[8] 31.22
## y1[9] 26.03
## y1[10] 21.45
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.2 seconds.
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 finished in 0.2 seconds.
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 finished in 0.2 seconds.
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 finished in 0.2 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 0.3 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -0.61 -0.24 1.66 1.44 -3.63 1.27 1.01 363 311
## b0 104.90 104.81 1.81 1.73 102.10 108.02 1.01 1072 800
## b1 0.03 0.03 0.00 0.00 0.03 0.03 1.01 1148 972
## b2 0.19 0.19 0.00 0.00 0.19 0.20 1.01 1056 868
## s 0.60 0.58 0.12 0.10 0.46 0.80 1.03 203 199
## m0[1] 57.76 57.76 0.22 0.21 57.41 58.12 1.00 1434 1357
## m0[2] 60.86 60.86 0.21 0.21 60.53 61.19 1.00 1762 1443
## m0[3] 38.52 38.50 0.37 0.35 37.94 39.12 1.01 1226 1296
## m0[4] 60.80 60.80 0.21 0.21 60.47 61.13 1.00 1745 1443
## m0[5] 32.72 32.72 0.20 0.19 32.40 33.05 1.01 1845 1305
## m0[6] 46.81 46.81 0.19 0.18 46.49 47.13 1.00 1623 1217
## m0[7] 61.41 61.40 0.21 0.21 61.06 61.75 1.00 1959 1478
## m0[8] 61.41 61.41 0.21 0.21 61.07 61.75 1.00 1952 1490
## m0[9] 22.93 22.93 0.25 0.24 22.52 23.33 1.01 1457 1228
## m0[10] 28.87 28.85 0.32 0.31 28.36 29.39 1.01 1185 1221
## m0[11] 37.07 37.07 0.19 0.17 36.76 37.38 1.00 1932 1245
## m0[12] 39.17 39.17 0.18 0.17 38.86 39.48 1.00 1911 1251
## m0[13] 24.71 24.71 0.24 0.24 24.31 25.10 1.01 1512 1275
## m0[14] 49.48 49.47 0.20 0.19 49.15 49.82 1.00 1509 1161
## m0[15] 23.98 23.98 0.25 0.24 23.58 24.37 1.01 1488 1279
## m0[16] 57.04 57.03 0.30 0.28 56.56 57.54 1.00 1634 1237
## m0[17] 54.01 54.01 0.22 0.21 53.66 54.37 1.00 1397 1266
## m0[18] 21.96 21.97 0.26 0.24 21.55 22.37 1.01 1431 1227
## m0[19] 44.47 44.47 0.19 0.17 44.15 44.77 1.00 1733 1317
## m0[20] 60.73 60.73 0.21 0.21 60.39 61.06 1.00 1731 1431
## y0[1] 57.77 57.79 0.65 0.61 56.71 58.84 1.00 1883 1838
## y0[2] 60.86 60.84 0.62 0.58 59.86 61.87 1.00 1791 1629
## y0[3] 38.53 38.51 0.72 0.67 37.39 39.69 1.01 1756 1540
## y0[4] 60.80 60.80 0.63 0.60 59.76 61.80 1.00 1816 1810
## y0[5] 32.75 32.73 0.64 0.62 31.70 33.83 1.00 2045 1774
## y0[6] 46.81 46.80 0.65 0.62 45.73 47.83 1.00 1818 1847
## y0[7] 61.41 61.40 0.64 0.61 60.37 62.51 1.00 2014 1434
## y0[8] 61.38 61.37 0.65 0.61 60.34 62.45 1.00 2002 1780
## y0[9] 22.91 22.90 0.68 0.61 21.83 24.02 1.00 1894 1497
## y0[10] 28.86 28.85 0.71 0.67 27.69 30.01 1.00 1545 1772
## y0[11] 37.06 37.06 0.64 0.62 36.00 38.05 1.00 2076 1829
## y0[12] 39.16 39.16 0.65 0.62 38.11 40.19 1.00 1858 1756
## y0[13] 24.68 24.69 0.67 0.62 23.56 25.78 1.00 2124 1587
## y0[14] 49.47 49.45 0.64 0.61 48.45 50.52 1.00 1831 1972
## y0[15] 23.99 23.98 0.66 0.63 22.94 25.07 1.00 1698 1654
## y0[16] 57.04 57.04 0.69 0.66 55.94 58.16 1.00 1964 1970
## y0[17] 54.02 54.01 0.64 0.64 52.99 55.10 1.00 1978 1902
## y0[18] 21.97 21.97 0.69 0.66 20.85 23.06 1.00 1853 1700
## y0[19] 44.47 44.47 0.65 0.62 43.41 45.54 1.00 2133 1659
## y0[20] 60.72 60.70 0.66 0.64 59.66 61.80 1.00 1942 1701
## m1[1] 28.87 28.85 0.32 0.31 28.36 29.39 1.01 1185 1221
## m1[2] 57.76 57.74 0.29 0.27 57.29 58.23 1.00 1675 1279
## m1[3] 61.06 61.06 0.21 0.20 60.73 61.39 1.00 1802 1452
## m1[4] 56.15 56.14 0.22 0.21 55.79 56.51 1.00 1400 1299
## m1[5] 49.17 49.17 0.20 0.19 48.84 49.50 1.00 1521 1161
## m1[6] 42.24 42.23 0.18 0.17 41.93 42.54 1.00 1821 1312
## m1[7] 35.99 35.99 0.19 0.18 35.68 36.30 1.00 1923 1187
## m1[8] 30.56 30.56 0.21 0.21 30.23 30.91 1.01 1772 1362
## m1[9] 25.92 25.92 0.24 0.23 25.54 26.30 1.01 1559 1279
## m1[10] 21.96 21.97 0.26 0.24 21.55 22.37 1.01 1431 1227
## y1[1] 28.86 28.85 0.71 0.68 27.74 29.99 1.00 1794 1529
## y1[2] 57.77 57.77 0.66 0.64 56.71 58.84 1.00 1984 1519
## y1[3] 61.07 61.05 0.64 0.61 59.98 62.11 1.00 2060 1563
## y1[4] 56.14 56.13 0.66 0.62 55.07 57.21 1.00 1919 1341
## y1[5] 49.18 49.21 0.65 0.62 48.09 50.20 1.00 2206 1883
## y1[6] 42.23 42.23 0.64 0.59 41.16 43.27 1.00 1941 1507
## y1[7] 36.03 36.02 0.63 0.61 35.00 37.06 1.00 1995 1914
## y1[8] 30.57 30.56 0.65 0.62 29.52 31.58 1.01 2077 1832
## y1[9] 25.92 25.94 0.65 0.65 24.86 26.97 1.00 1988 1584
## y1[10] 21.97 21.97 0.66 0.63 20.88 23.07 1.00 2006 1938
